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Abstract. We report a new approach for matching regular surfaces in a Rie- 
mannian setting. In this setting a deformation vector field is defined on the 
surface directly and the motion of the surface does not involve motion of the 
ambient space. This is in contrast to LDDMM, in which a prescribed defor- 
mation of the ambient space drags along an embedded surface. In particular, 
this approach performs matching using an inner metric defined on the space of 
surfaces. We also show how to discretize the corresponding geodesic equation 
and compute the gradient of the cost functional using finite elements. 



1. Introduction 

We introduce a new paradigm for image registration in computational anatomy, 
in which the image characteristics of regular surfaces embedded in R 3 are trans- 
formed directly from a template to a target image, without requiring any deforma- 
tion of the ambient space. 

The field of computational anatomy concerns itself with the study and classi- 
fication of the variability of biological shapes, including their statistical variance. 
Since the space of all shapes is inherently nonlinear, the usual methods of linear 
statistics cannot be applied. In particular, the addition of two surfaces cannot be 
meaningfully defined. One way to overcome this difficulty is to introduce a Rie- 
mannian structure on the space of shapes, which locally linearizes the space and 
allows the development of statistical methods that are analogous to the linear case. 
This approach was taken, e.g., in [6]. In the Ricmannian setting, the average of two 
shapes may be defined as the middle point of a geodesic joining these two shapes. 
In a similar way one may define the corresponding geodesic mean of a collection of 
n shapes. 

One class of shapes which are of interest in computational anatomy consists 
the surfaces embedded in R 3 . The cortical surface, the surfaces of hippocampi, 
thalami, putamen and caudate are all examples of shapes which are represented as 
two dimensional surfaces in R 3 . This is also an example, where the Riemannian 
setting may be applied to study collections of shapes. 

The currently prevalent method for comparing anatomical shapes in the Rie- 
mannian setting is the method of large deformation diffcomorphic metric matching 
(LDDMM), based on the deformable template paradigm of Grenander [8]. In this 
setting, a template shape is matched to a target shape by finding a transformation in 
a suitable group of deformations of the ambient space that transforms the template 
into the target. This approach has been systematically developed in [3, 4, 7, 14, 15] 
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and applied to various problems in computational anatomy. Registering two sur- 
faces in this framework involves finding a diffcomorphism of the whole space, which 
transforms one surface to the other. 

In this paper we propose a different way of defining a Riemannian structure on 
the space of surfaces, which also provides the full range of tools for nonlinear sta- 
tistics. Unlike LDDMM, our approach applies the transformation theory of surface 
registration directly to the image, without any deformation of the ambient space. 
Our approach to characterizing a deformation is intrinsic to the surface, rather 
than resulting from a transformation of the surface induced by a deformation of 
the ambient space in which the image is embedded. For this reason, the Riemann- 
ian metrics used here are called inner metrics as opposed to the outer metrics used 
in LDDMM, where deformations are imposed via the ambient space. 

Inner metrics on planar curves were introduced and studied previously in [11, 12, 
16]. Recently several examples of inner metrics were studied on surfaces and higher 
dimensional hypersurfaces in Euclidean space in [1, 2]. The numerical implementa- 
tion of inner metrics for matching differs from LDDMM, because the metric on the 
tangent space at each surface depends nonlinearly on the surface. This means the 
metric will change adaptively as one moves around in shape space. This adaptive 
property of the inner metrics is in marked contrast to LDDMM, where the outer 
metric is defined on the diffeomorphisms and projected down to the shapes, so it 
doesn't depend on the particular shape. 

The outline of the paper is as follows. In Sec. 2 we review the registration prob- 
lem for surfaces and recall how it is solved using outer metrics in LDDMM. Then 
we present the approach via inner metrics and point out the differences between the 
two methods. For definiteness, we will concentrate our attention on the Sobolev 
metric of order one. In Sec. 3 we discuss how to discretize and implement the 
geodesic equations for this metric and how to solve the registration problem via 
geodesic shooting. Finally, in Sec. 4 we show how this metric performs in some 
examples using synthetic data. 

2. The Mathematical Formulation 

We are dealing with the registration of parametrized regular surfaces. Such a 
surface is given by a smooth function q : M — > M 3 from a model surface M into 
the Euclidean space M 3 . We will consider different choices of the model surface M 
in this paper: the plane sheet M = [0, 1] x [0, 1], the cylinder M — S 1 x [0, 1] and 
the torus M = S 1 x S 1 . Another interesting choice would be the sphere M = S 2 , 
which however is not considered in this paper. Treatment of the sphere is more 
challenging, because it cannot be represented by a single global coordinate chart. 
We require the parametrization of the surface q to be regular in the following sense: 
at each point x € M the partial derivatives g^-, are required to be linearly 
independent. We will denote the space of all such surfaces by Y . 

2.1. Registration with LDDMM. In the LDDMM framework, the registration 
of a template surface q to a target surface cftarg involves finding a curve ip t of 
diffeomorphisms of the ambient space ]R 3 , whose deformation carries the template 
surface to the target surface. Mathematically, one constructs these deformations 
using time dependent vector fields v t {y) 1 which generate ip t as their flow, i.e. 

(1) d t (p t = v t (tp t ) . 
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Figure 1 . For inner metrics the vector field governing the defor- 
mation is defined directly on the surface (left picture) . In contrast 
the LDDMM approach defines a deformation vector field on all 
of M 3 (right picture) . The latter vector field deforms the ambient 
space and in the process induces a deformation of the surface. 



The registration problem consists of finding a vector field, which minimizes the 
following sum of a kinetic energy and a matching term 

(2) E{v t ) = ^J^ |H|yttt + ^dfaiGfo),«targ) • 

The kinetic energy is usually measured using a norm \\.\\v defined on a reproducing 
kernel Hilbert space V of vector fields on M 3 with kernel K . The norm is then given 
by \\ u \\v — I u ' K~ l *udx. We will discuss possible choices of the matching term 
in Section 2.5. 

It is possible to reduce the complexity of the problem, since one can show that 
minimizing vector fields v t have to obey an evolution equation. This enables us to 
describe the whole vector field vt by knowing only its value vo at time t = 0. The 
equations 



(3a) d t q t (x) = v t (q t (x)) 

(3b) d t p t (x) = (Dv t {q t (x))) T Pt (x) 

(3c) v t (y)= p t (x)K(y-q t (x))dx 

JM 



are given in terms of a momentum p t , which lives on the surface. The momentum is 
convolved with the kernel in (3c) to reconstruct the minimizing vector field, which 
defines the deformation of the ambient and drags the surface along in (3a) . Details 
about LDDMM and the reformulation as evolution equations can be found in [3, 17]. 

2.2. Registration with Inner Metrics. We propose to use a different approach, 
described from a mathematical point of view in [1, 2]. In this approach we describe 
the deformation of the surface directly, without assuming an underlying deformation 
of the whole space. In our approach we will replace (3a) by 

(4) d t q t {x) = u t (x) , 

where u t (x) g C°°(M, M 3 ) is a time dependent vector field, defined only on the 
surface. Note the difference between the vector field v t (y) which is defined on R 3 
and u t (x), which is defined on the model space M, c.f. Fig. 1. 
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In the new framework the registration problem still consists of minimizing the 
following sum of a kinetic energy term and a matching functional 



The kinetic energy is measured via an inner product (., .) qt on the space of vector 
fields along the surface. This inner product can and usually will depend nonlincarly 
on the surface qt itself. This is another difference with the LDDMM framework, 
where the kinetic energy of the vector fields didn't depend on the surfaces that were 
matched. There is a whole variety of inner metrics, one can choose from. We will 
concentrate in this paper on an i? 1 -type metric, which will be discussed in Sec. 2. 4. 
Again we postpone the discussion of the matching term to Sec. 2.5. 

2.3. Geometry. Both the LDDMM approach and the inner metrics can be seen as 
a special case of a constructions in the general framework of Riemannian geometry. 
In the case of inner metrics the collection of inner products (., .) q defines a Riemann- 
ian metric on the space of all surfaces. In Riemannian geometry, curves, which 
minimize the energy J Q (ut,ut)q t dt for fixed endpoints qa,qi are called geodesies. 
We can see that minimizers of (5) have to be geodesies in the space 5?. 

In LDDMM the inner product on the space of vector fields V also defines a 
Riemannian metric, this time on the group of diffeomorphisms. Therefore the 
minima of the registration problem (2) generates geodesies in the diffeomorphism 
group. 

Why is it advantageous to work in a Riemannian setting? In this setting, the 
minima of the matching energy are geodesies, so one may describe the nonlinear 
space of shapes in terms of the initial velocity or momentum of a geodesic, which 
is an element in a linear vector space. This is possible, because geodesies obey an 
evolution equation like (3). Using the initial velocity, which encodes the whole 
solution of the registration problem, we are able to view the space of surfaces from 
the template surface qo as a linear space. This enables us to use statistics, compute 
average surfaces and measure distances. 

2.4. H 1 -type Metric on Surfaces. We will consider a metric on the space of 
surfaces, which is the analogue of the i/^-norm ||/||^i = / R2 |/0c)| 2 + \Vf(x)\ 2 dx 
for functions on M 2 . We will replace functions on R 2 by vector fields living on the 
curved surface q and adjust the definition of the if -norm to take into account the 
curved nature of the surface q. Since q is a surface in R 3 , we can measure angles 
and distances of vectors tangent to the surface, using the Euclidean inner product 
on M 3 . At each point q(x) of the surface we also have a canonical basis for the plane 
tangent to q, given by the vectors J^, J^. We denote the inner product induced 
on the tangent plane to the surface by g. This inner product has the following 
coordinate matrix with respect to the basis 4^ , A ■ 



(5) 





-l 




We denote by (<7 4J ) the inverse matrix of (gij) and by vol(g) the volume density of 
the surface q with respect to the metric g. 
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In this paper we will use the H -type inner metric on surfaces defined by the 
expression 

3 

(7) (u, v) q : = y2 / u k v k + a 2 g (grad 9 (u fc ), grad 9 (v fc )) vol(g) dx 

One reason to use this generalization of the if ^metric is, that this metric is invari- 
ant under reparametrizations of the surface and only depends on the image q(M) 
as a subset of M 3 , in a similar way as the length of a curve in two dimensions only 
depends on the image of the curve, and not on a particular parametrization. This 
is necessary, if one wants to match unparametrized surfaces. For this task it is 
possible to use the same framework with this metric, only the matching term has 
to be chosen to beinvariant under reparametrizations. 

The constant a, which appears in the metric, is a parameter, which has to be 
chosen for each problem. It represents the characteristic length scale, at which 
deformations take place. Another interpretation of a is the scale across which 
the momentum is smoothed, when passing from momenta to velocities. It can be 
compared with the kernel size in LDDMM. 

Other choices for the metric are possible. One could use Sobolev type metrics 
of higher order as in [2] or multiply the components of the metric with a func- 
tion depending on geometric quantities of the surface, like area, mean or Gaussian 
curvature as was done in [1]. 

2.5. The Matching Functional. There are different possible choices for the match- 
ing term. In this paper we will use the squared L 2 -distance 

(9) d(q ,qi)= / \q (x) - qi{x)\ 2 dx. 

Jm 

Since we are dealing with parametrized surfaces this is a natural choice for the 
matching functional. 

When matching unparametrized surfaces, natural choices of the matching func- 
tional would include currents, see [7], or one could use the reparametrization frame- 
work of [5] . 

3. Discretization 

In this section we will describe how to discretize the variational problem (5) and 
compute the optimal path between two surfaces. Starting with an initial guess for 
the velocity Uq, we will use a gradient descent scheme 

to converge towards the initial velocity of the optimal geodesic. The discretization 
thus consists of two parts: 

• compute the geodesic, given the initial velocity to evaluate E(uq) 

• compute the gradient V Ua E(u l ) to update the initial velocity. 

We show how to discretize the geodesic equation in Sec. 3.1 and how to compute 
the gradient in Sec. 3.2. 
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3.1. The Geodesic Equation. We discretize the time-evolution of the surface 
q(t) using the explicit Euler method 

(10) q i+1 =qt + AtUi , 

where qi — q(iAt) is the discretized version of the curve and At = 1/N is the 
time step, if we divide the interval [0, 1] into N parts. To compute Ui we note 
that a geodesic is a critical point of the energy E(ui) = ^ Si=o ( u » i u i)qi> 1 - c - 
\7 Ui E{ui) = 0. Following [9] we introduce the Lagrangian multiplier pi in the 
discrete variational principle 

(11) E(u ,...,u N -!) = ^ —(ui,Ui) qi + (pi,q i+ i - qi- AtUi) L 2 . 

i=0 

and take variations. From variations in Uj we see that Pi is the momentum dual to 
the velocity Ui in the sense that (u,, Sui) qi = (pi, Siti) L 2 and we obtain the evolution 
equation for Mj in the form 

/ 5£ 

(12) (u i+1 ,5q i+1 ) qi+1 = (ui,6q i+1 ) qi + At ^ — (u i+1 ,u i+1 ;qi +1 ),8q i+1 ^ 
with 5qi+i an arbitrary variation. Here we use the notation 

(13) i(u,v;q) = ^{u,v) q . 

We denote by || the variational derivative of £(u, v; q) with respect to the variable 
q, defined via 

/ 51 \ £(u,v;q + h5q)-£(u,v;q)) 

(14) {—(u,v;q),Sq) = lim . 

\oq I L 2 fr->o h 

Equation (12) is an implicit time step for u i+ i, since u i+ \ appears on the right 
hand side in a quadratic term. To make computations easier and avoid having to 
solve a nonlinear equation, we changed the right hand side to an explicit Euler time 
step 

(15) (u i+ i,6q) qi+1 = (ui,6q) qi + At(^-{ui,Ui,qi),5q^ 

3.2. Computing the Gradient. Given (10), (15) for the evolution of a geodesic 
we again use method of adjoint equations from [9] to compute the gradient of the 
energy with respect to the initial velocity. The resulting equations for the variables 
Ui,Vi have to be integrated backwards in time 

{ui,5ui) q% = (u i+ i,8ui) qi + At(v i+ i,Sui) qi+1 +2At( — (u i ,5u i ;q l ),u i+ i) 

I 81 

(16) (vi,5 qi ) qz = (v i+1 ,6qi) qi+1 + 2(—(u i+1 -u i ,u i ;q i ),Sq i 

I 5 2 £ 

with the initial conditions 

(17) u N = (v N ,5q N ) qN = -^{Qn - tftarg'^iv) 
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Figure 2 . Samples are shown from a geodesic in the space of sur- 
faces between a straight cylinder and a bent cylinder with ripples, 
at time points t = 0,0.3,0.5,0.8,1. The color encodes the Eu- 
clidean length of the deformation vector field at each point of the 
surface. 

at time t = 1. The gradient is then given by 

(18) V Uo E(u ) =uq-u ■ 



We implemented the geodesic and adjoint equations (10), (15) and (16) in Python 
using the finite element library FEniCS [10]. All model manifolds ([0,1] x [0,1], 
S 1 x [0,1], S 1 x S 1 ) were modelled on the rectangle [0,1] x [0,1] with periodic 
boundary conditions prescribed where neccessary. The domain was subdivided into 
a regular triangular mesh, on which Lagrangian finite elements of order 1 were 
defined. 

In the first example we apply our method to compute the geodesic path between 
two shapes, which includes both large and small deformations. The template shape 
is a straight cylinder of height 1 and radius 0.25, which is discretized using a regular 
triangular mesh of 2 x 30 x 30 elements. The target shape is a cylinder, which is 
bent by 90° and has 5 small ripples added to it along the vertical axis. Compared 
to the bending the ripples constitute a small and local deformation of the shape. 
The target shape is discretized in the same way as the template. We use a = 0.6 as 
the length scale parameter and 10 time steps for the time integration. The gradient 
descent takes 80 steps to converge to an L 2 -error of 0.008. We can see in Fig. 2 
that both the large and the small deformations are captured by the geodesic. 

In the second example we want to illustrate the curved nature of shape space. 
To do so we pick three asymmetric tori, lying in different positions in space. Each 
two tori differ by a composition of two rigid rotations. We compute the geodesies 
between each pair of tori to measure the angles and side lengths of the triangle 
with the tori as vertices and the geodesies as edges. By comparing the sum of the 
angles with n one can estimate, whether the curvature of shape space along the 
plane containing the triangle is positive or negative. We measured a = 33.766° , 
P = 34.802° and 7 = 34.675° . The sum a + /3 + 7 = 103.243° is smaller than 
180° , which indicates that the space is negatively curved in this area (c.f. [12, 
Sec. 5.4]). In negatively curved spaces geodesies tend to be attracted towards a 



4. Numerical Experiments 
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Figure 3. This figure shows a geodesic triangle in the space of 
surfaces with asymmetric tori as vertices. The tori along the edges 
are the middle points of the geodesies connecting the vertices. One 
can see that the shapes tend to shrink along the geodesies before 
expanding again towards the ends. This effect implies negative 
curvature in this region of shape space. 

common point. In this example the geodesies are attracted towards the surface, 
which is degenerated to a point. We can see in Fig. 3 that the midpoints of the 
geodesies between the vertices are slightly shrunk. This is another indication for 
the negatively curved nature of the space. 

In the third example we show, that our framework is indeed capable to do non- 
linear statistics on shape space. We generate five sample shapes and compute the 
mean shape between them. The five shapes are cylindrical vases with an open top 
and bottom, discretized again using a triangular mesh of 2 x 30 x 30 elements. As 
the initial guess for the mean we use a straight cylinder. First we register this initial 
shape to the five target shapes, compute the average of the initial velocities and 
then shoot with this average velocity to obtain a next guess for the mean shape. 
We iterate this procedure until the average velocity is close to zero. This method of 
computing the Karcher mean was proposed in [6] . After four iterations we obtained 
an average velocity with norm 0.006. As can be seen in Fig. 4 the average shape 
indeed combines the characteristics of the five shapes. 

5. Conclusions 

In this paper we propose a new approach to match regular surfaces in a Rie- 
mannian setting. Although this approach has been studied from a mathematical 
point of view in several papers, including [1, 2, 13], so far it hasn't been applied to 
problems in computational anatomy. The aim of this work is to argue, that this is 
a promising approach, which is worthwhile to be studied further. 

In contrast to LDDMM, where surfaces are deformed via a deformation of the 
ambient space, in this approach the deformation is prescribed directly on the sur- 
face, while the ambient space stays constant. Because of this we call this approach 
matching with inner metrics as opposed to LDDMM, which can be described as 
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Figure 4. In this figure we show the Karcher mean of five vase- 
shaped objects. The mean shape, which is displayed in the center 
of the figure is computed using an iterated shooting method. The 
colored regions on the averaged shapes encode the Euclidean length 
of the initial velocity of the geodesic, which connects each shape 
to the mean. The color of the mean was chosen for artistic purposes 
only. 

matching with outer metrics. We show how to discretize the geodesic equation and 
how to compute the gradient of the matching functional with respect to the initial 
velocity. In the last part of the paper we present numerical results on synthetic 
data of different topologies, which demonstrate the versatility and applicability of 
our method. 

At the present we applied this metric to match parametrized surfaces, which is 
an unwelcome restriction in applications. This is not a restriction of the framework 
itself, but only of the matching functional. By choosing a matching functional, 
which is independent of the parametrization of the surface, one can apply the same 
framework to unparamctrized surfaces. In future work we plan to implement this 
capability and test the method on real anatomical data. 
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